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ABSTRACT 

An analytical approximation to periodic orbits in the circular restricted three- 
body problem is provided. The formulation given in this work is based in calcu- 
lations known from classical mechanics, but with the addition of the necessary 
terms to give a fairly good approximation that we compare with simulations, re- 
sulting in a simple set of analytical expressions that solve periodic orbits on discs 
of binary systems without the need of solving the motion equations by numerical 
integrations. 
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1. INTRODUCTION 

The majority of low-mass main-sequence stars seem to be grouped in multiple systems 
preferentially binary (Duquennoy & Mayor 1991, Fischer & Marcy 1992). These systems 
have attracted more attention since the discovery that many T-Tauri and other pre-main- 
sequence binary stars possess both circumstellar and circumbinary discs from observations 
of excess radiation at infrared to millimeter waves and direct images in radio (Rodriguez 
et al.l998; for a review see Mathieu 1994, 2000). The improvement in the observational 
techniques allows, on the other hand, the testing of theories on these objects. Moreover, 
extrasolar planets have been found to orbit stars with a stellar companion, e.g., 16 Cygni B, 
T Bootis, and 55 p Cancri (Butler et al.l997; Cochran et al. 1997). All these facts make the 
study of stellar discs in binary systems, as well as the possibility of stable orbits on them 
that could be populated by gas or particles, a key element for better understanding stellar 
and planetary formation. 
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In particular the study of simple periodic orbits of a test particle in the restricted three- 
body problem results a very good approximation to the streamlines in an accretion disc in 
a very small pressure and viscosity regime, or for proto-planetary systems and its debris 
(Kuiper belt objects). Extensive theoretical work has been done in this direction (Lubow & 
Shu 1975; Paczyhski 1977; Papaloizou & Pringle 1977; Rudak & Paczyhski 1981; Bonnell & 
Bastien 1992, Bate 1997; Bate & BonneU 1997). 

Several studies arc directed to find the most simple and important geometric character- 
istic of the discs which is the size of both circumstcllar and circumbinary discs, going from 
analytical approximations: Eggleton (1983), who provides a simple analytical approximation 
to the Roche lobes; Holman & Wiegert (1999) radii in planetary discs; Pichardo, Sparke & 
Aguilar (2005, PSA05 hereafter) radii in eccentric binary stars. 

Based in the approximation studied in classical mechanics to periodic orbits given by 
Moulton (1963), we provide in this work a fairly good approximation by calculating the 
necessary terms, for the case of circular binaries, to periodic orbits. This formulation gives 
not only the radii of the circumstcllar and circumbinary discs, but gives a good description 
of any periodic orbit at a given radius for any of the discs, in the form of an analytical 
approximation. In this manner, the set of equations we provide can be used to find any 
periodic orbit or its initial conditions to run it directly in simulations of the three-body 
restricted problem, or as a good approximation to initial conditions in the eccentric case as 
well, from the point of view of particle or hydrodynamical simulations, without the need of 
solving the restricted three body problem numerically. 

The low viscosity regime described by the periodic orbits representation, is important 
in astrophysics since it gives an idea of how bodies would respond only infiuenced by the 
effects of the potential exerted by the binaries, and they could permit to link the periodic 
orbits to physical unknown characteristics of the discs like viscosity that needs to be artifi- 
cially introduced to hydrodynamical codes. It also might permit the calculation of physical 
characteristics like dissipation rates, etc. Studies like these result easier given in the form of 
approximated analytical expressions, instead of the usual numerical approximations to the 
periodic orbits since an analytical formulation is much faster to add to any hydrodynamical 
or particle code that requires the characteristics of a given family of periodic orbits, or sim- 
ply their initial conditions in a more precise form than assuming Keplerian discs until the 
approximation fails. The periodic orbits show regions where the orbits are compressed (or 
decompressed), these regions could trigger density fiuctuations maybe able to drive material 
to form important agglomerations in the discs that, depending on their positions on the 
disc and the density, could give origin to planets. An analytical approximation with a given 
density law for the discs would allow these kind of studies to be much faster and easier. 
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In Section [2l we show the strategy and equations used to find the approximation to the 
periodic orbits in circumstellar discs. In Section [3l we present the equations approximating 
the periodic orbits in circumbinary discs, and an analysis of stability in Section I3.2[ A 
comparison of the application of the formulation and numerically calculated periodic orbits 
is given in Section HI Conclusions are presented in Section [51 

2. Methodology 

In general, orbits around binary systems are calculated by approximate methods. Nu- 
merical calculations are the most extended. In this work we are going to use other known 
method to search for an analytical approximation that will provide us of solutions for any 
orbit in circumbinary or circumstellar discs of the circular problem. The method is based 
on a perturbative analysis. In this case we choose some terms in the expanded equations of 
motion, which are relevant to the solution of the problem. The original equations are approx- 
imated with a set of equations that, in many cases, have an analytical solution. Comparing 
with the numerical approach, this one has the advantage that the expressions are completely 
analytical and simple. We compare our calculations with numerical studies in order to pick 
as many terms in the expansion as we need to obtain the solution. The method can be used 
as a way to find a quick answer and for studies in which an analytic expression makes the 
problem manageable. 

2.1. Periodic Orbits in Circumstellar Discs 

The analysis is restricted to the orbital plane, in polar coordinates the equations of 
motion are given by, 

r — rip^ = Fr, 
ri! + 2rip = F^, (1) 

where the origin of the system is located at the position of the star with mass Mi (hereafter 
we call this star 5*1 ). It is important noticing that 5*1 can represent either the primary or 
the secondary star. Here, r and are the radial and angular coordinates respectively, and 
the forces along this directions {Fr,F^) have the form 
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that represents the perturbing potential (divided by GM2) due to the star with mass M2 
(hereafter, star 5*2). $p could be also expressed in terms of the Legendre Polynomials (Murray 
& Dermott, 1999). Here R and ^ are the coordinates of S2 and r and ip locate the particle 
that is perturbed, P. 

The non-perturbed state corresponds to the case with $p = 0. Thus, the non-perturbed 
orbits are conies around 5*1, which are perturbed by the presence of the other mass. The 
set of equations ([1]) cannot be solved analytically, then we expand the perturbation given in 
equation ([2]) using r/i? as a small parameter. This parameter at some points in an orbit that 
begins around a resonance is large, so this approach is not valid. This expansion improves 
in precision the closer the particle is to Si. Thus, equations ([1]) can be written as. 
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(3) 



which is equivalent to the expansion (6.22) in Murray & Dermott (1999). 

In this analysis the star 5*2 orbits S\ describing a circular orbit. We consider the mass 
of the discs negligible compared with the mass of the binary system. Thus, the orbit of 5*2 
is given by 
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R = a, 

^ . y£(K+M(,_,„).0(*-*„), (4) 

where a represents the fixed distance between the stars, and fl is the angular velocity of 5*2 
from the third law of Kepler. At t = to the star is in the reference line, which for a circular 
orbit can be any radial line that begins on Si. The non-perturbed trajectory of the particle 
P, is a Keplerian orbit around 5*1, in this manner the polar coordinates are expressed. 



where a„p and uj have the same meaning of its counterpart in equation HJ but relating P and 
Si. 

The perturbed position of the particle P can be written, 

r = a„p(l + rp), (6) 

ip = uj(t- to) + tpp, (7) 

where the terms {anpTp) and ipp are the corrections in the trajectory due to the presence of 
5*2. We define now the parameters rUo and r as follows, 

mo 



UJ — Vt' 

T={uj-Q){t-to), (9) 

here rUo is a small parameter for r <^ R. This condition is required for the expansion of the 
equations of motion (equations [3]) to be useful, r is the angle that locates the particle P, 
measured from the line connecting the stars, for the non-perturbed orbit. The first step is 
to replace t for r in equation ([3]), using equation ([9]). Equations ([6] and [7]) are substituted in 
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equations ([3]) taking into account that and tpp depend on r. The resulting expressions are 
finally divided by anp(^)^- The mass is given in units of Mi + M2, and Mi + M2 = 1, then, 



■ \2 (1 + ^o) 

Tp - {1 + rp){l + mo + ^p) + (i^^^)2 



mlM2 



l + rp)[(l + 3cos2(r + V^p)) 



3 

+ -Mpmo(l + rp)(3cos(r + ijjp) + 5cos3(r + ^pp)) 



9 5 , .35 
8 + 2 

+...] , (10) 



+M^ml{l + rp)2(- + - cos 2(r + t/^p) + — cos4(r + ^p)) 



(1 + rp)'iljp + 2rp(l + mo + z/^p) 
-^l^(l + rp)[3sin2(r + ^p) 



3 

+ -Mpmo(l + rp)(sin(r + ipp) + 5sin3(r + ijjp)) 

5 7 
+-M2m2(l + rp)2(sin2(r + ^p) + - sin4(r + ^/-p)) 

+...] . (11) 

Equations (fTOl and [TT]) are dimensionless. They will be solved for Vp and ipp in terms 
of M2 (Ml = 1 — M2), r and rrio. The value for M2, characterizes the stellar masses of 
the system, r defines the angular position of the point in the original non-perturbed orbit, 
and nio depends on the radial position of the non-perturbed trajectory Onp, in units of the 
separation between the stars, a. An analytic solution is found if an expansion for Vp and ipp 
in powers of nio is made as follows, 

rp = J:jL2rp,{M2,T)mi, (12) 
iPp = ^%2MM2,r)mi. (13) 

Specifically, for periodic orbits, it is necessary to fulfill the next conditions, 

rp,{r + 27i) = rp^.(r), 

^p^(r + 27r) = iPp^ir). (14) 
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An useful assumption is that the angular position of the orbit, from the line that connects 
the two stars, is not perturbed by 5*2, a fact necessarily true, which can be demonstrated 
using symmetry arguments. We restrict the perturbed orbits to be symmetric with respect 
to the line joining the stars by. 



V'p,(0) = 0, (15) 

this condition is general enough since there are no restrictions provided by observations 
about the orientation of the symmetry axis in the plane of the discs with respect to the line 
joining the stars. 

Equations ([12] and [13]) are substituted in equations ([10] and [11]) and expanded in powers 
of rUo, ending with a polynomial in rrio equal to zero. The only solution for such an equation 
is that every coefficient is zero. Each coefficient associated to a given power, is a second 
order differential equation. The solutions of the equations for with k = 2,3, 4, are solved 
imposing the conditions given in equations (HM and ]T5l) and are given by 
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(18) 
(19) 
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with 
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256 



MlH + 



1215 
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(22) 



It is worth noticing that for each order (in m^) we obtain a system of two differential 
equations that produce two integration constants that can be calculated by solving the set 
of equations of the next higher order, by imposing the conditions in equations [H] and [151 

The parameter H comes from the approximation a„p = rrioH, where we will consider H 
constant, justified by the fact that H depends weakly on nio (Moulton 1963). 

Using equations f|T6l [T71 [T8l [T9| [20] and [21]), plus equations f[T2] and [13]) in equations ([6] 
and [7]) will let us find the perturbed trajectory of the particle P. Since a„p is proportional 
to rrio, the radial corrections are of order and the angular corrections are of order m^. 
Thus, equations ([6] and [7]) represent coordinates depending on time, which means velocity 
and accelerations can be obtained by direct differentiation. It is worth noticing that the 
orbits are given in the reference system where the stars are fixed in the X axis. An example 
of the orbits and rotation curves produced by this model are given in Section [H 



To calculate the total radial size of discs we will assume we are in the regime of low 
pressure gas. The permitted orbits for gas to settle down there, are all those that do not 
intersect themselves or any other (Paczyhski 1977, PSA05). In this paper we will use the same 
approximation to find discs radii except that the intersections are found analytically (other 
analytical approximation to the radial extension for circumstellar discs for any eccentricity 
can be found in PSA05). First, we choose a given M2, to fix the binary system. A physical 
intersection occurs when 



where is the radius of the orbit given by equation and ipi spans the 27r-angular-range 
of the same orbit. The subindex i takes the values i = 1,2 for inner and outer consecutive 
orbits, respectively. We now look for a value for a„p and r that simultaneously satisfy 
equations f[23] and \2M . Note that equation f[23]) can be expanded in a series of Taylor in the 
variable a^p. Due to the assumption that we are interested in infinitesimally close orbits, a 
linear approximation is good enough. In this manner, equation f[25]) can be transformed to 
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The Radial Extension of Circumstellar Discs 



^2 - V^i = 



r2 - ri = 0, 



(23) 
(24) 



which is a function of M2,a„p and r. Also, using equation ([7j), the equation can be 
written as 



From equations ( fTTl [T9] and [2T1) we directly identify that all the terms are proportional to 
sin(A;r), which results proportional to sinr. Thus, sinr can be factorized in equation (1261) 
giving the solutions r = and r = vr. Equation fl26|) has other solutions but none allows to 
find a root for equation fl25l) . Taking (r = 0,r = n) and using any of these options in equation 
( l25l) . we end up with an equation for a„p. The solution of equation ( l25l) using r = 0, vr gives 
two values of a„p, which correspond to the intersections. We take the intersections at the 
smallest radius (that represent the maximum radii where gas particles would be able to settle 
down), thus the second intersection at larger radius are not useful in this analysis. In both 
cases the intersecting orbits are tangent to each other, and contiguous orbits between them 
intersect at an angle. Because we do not consider this kind of orbits as part of the disc, the 
disc naturally ends up in the orbits with smaller a„p. 

The minimum value for a^p comes from r = tt, this is, the intersection appears at 
the opposite side of the star 5*2. In Table [H the solution of equation fl25|) for different 
values of Mi, with r = vr, is given in the column named anp^r=n- The following two columns 
(r(r = 0),r(r = vr)) give the radial positions of the innermost intersecting orbit. The average 
of these values is given in the column < r >, the next column gives the same average but 
from PSA05, and the last column is the approximation to the Roche Lobe radius by Eggleton 
(1983) given by, 



^p,2 - ^p,l = 0. 



(26) 




0A9qt 



where 



(27) 



a 



a 



0.6qr + ln{l + q. 



1/3 



Qi = 



I-M2 

M2 
M2 



(28) 



(29) 
(30) 



q2 = 



1-M2 ' 



where i refers to either star Si with i = 1,2. 
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From Table [H < r > and < r >psA (for particular masses of the stars, Mi and 
M2 = 1 — Ml) coincides within a 20% of error for Mi < 0.5, decreasing to less than 10% for 
Ml > 0.5. As expected the last non-intersecting orbit is contained inside the Roche lobe, as 
we can see in the Figured] for M2 = 0.2. A good approximation for the Roche lobe is given 
using the approximation of Eggleton (1983). 



3. Orbits in Circumbinary Discs 



The method described in Section [2711 is followed in this case with few differences. First 
of all, the origin is now on the center of mass of the system. The perturbing potential 
takes the form, 



GMi 



1-^^ cos('?/'-^/'i) + (— )^ 



-1/2 GM2 
+ 



l + ^cos(^-^i) + (:^)^' 



-1/2 



(31) 



where Ri = M2R, R2 = MiR are the distance to the star Si and S2 from the origin, 
respectively. It is important to mention that 5*2 is located to the left of the center of mass, 
and that it can represent either the secondary or the primary star, a is the distance between 
the stars, r and ip — tpi are the coordinates for the particle P that trace the circumbinary 
orbit, where the angle is measured from the radial vector which aims to Si. In this case, 
equation ([1]) can be used with the new definitions as follows. 



dr 



(32) 



Table 1: Sizes of the circumstellar discs. 



Ml 




r(r = 0) 


r(r = tt) 


< r > 


< r >PSA 


Rrl 


0.1 


0.1628 


0.1641 


0.1326 


0.1484 


0.125 


0.213 


0.2 


0.2077 


0.2113 


0.1690 


0.1901 


0.162 


0.268 


0.3 


0.2445 


0.2486 


0.1989 


0.2237 


0.195 


0.308 


0.4 


0.2795 


0.2825 


0.2276 


0.2550 


0.228 


0.344 


0.5 


0.3155 


0.3155 


0.2573 


0.2864 


0.257 


0.379 


0.6 


0.3543 


0.3489 


0.2901 


0.3195 


0.317 


0.414 


0.7 


0.3982 


0.3841 


0.3285 


0.3563 


0.350 


0.454 


0.8 


0.4507 


0.4236 


0.3766 


0.4001 


0.387 


0.501 


0.9 


0.5213 


0.4774 


0.4466 


0.4620 


0.426 


0.570 



(33) 



Thus, an expansion of equation ([T]) for ^ <^ 1 and <^ 1 can be developed, 



f — rij)^ + 



G(Mi + Ma) 




V^i)] + •••}, 



rip + 2rilj 




V^i)] + ... 



}, (34) 



where Mp = M1M2 = M2(l - M2). 

The third term on the left side of the first equation represents the contribution to the 
force, in the case that all the stellar mass is concentrated in the origin. The right side of 
both equations are the perturbing force due to the mass distribution between the stellar 
components. If the right side of equation (l34ll is equal to zero, we find the non-perturbed 
orbit in the form. 



Note that the definition of uj differs from the circumstellar case in the mathematical 
expression, however the meaning is the same, non-perturbed angular velocity of the particle 
P. Again, the perturbed position of P is given by equations ([H]and[7]). Equation changes 
to 




r 



a. 



(35) 



(36) 



T = 



{n-uj){t-to) 



(37) 



where Q = ^/^^^^ i 



is the angular velocity of the binary system, then r represents the 



angle between 5*1 and P. 



An useful parameter is given by 
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where /z is the parameter used to expand equation then /i is small as we are only taking 
into account circumbinary orbits far away from the stars. Changing t ^ t (equation [371) 
and expanding equation fl34p in a power series of /x^^^, we can write the equations of motion 

as 



rp-(l + rp)(-^ + V^p)' + 



1 



/i 



(2M2 - 



10/3 3 

|-[l + 3cos2(r + 7/'p)] + 



2(1 + r„ 



[3 cos(r + ifjp) + 5 cos 3(r + ipp 



(39) 



(1 + rp)?/Jp + 2rp(— Vtpp) 

1 — u 



Mp /iio/3 ^3 



3 (2M2 - l)/it. . 



|-sin2(r + 7/^p) + 



(1 + Tp) 



sin(r + ^p) + 5sin3(r + V^p)] + ...|. (40) 



In the same way as in Section [2711 differential equations are extracted at different orders 
if we expand rp and ipp as follows, 

rp = S-4rp,(M2,r)/x^/3, (41) 
^P = S~4V'p.(M2,r)/x^/l (42) 

The perturbed trajectories must be closed, then equations ([T^ are imposed. Also, the 
fact that the orbit is symmetric with respect to the line joining the two stars, allows to find 
the following condition. 



rp^(r = 0) = 0. 



(43) 
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Substitution of equations (jH] and |12]) in equations ([39] and HQ]) gives a couple of poly- 
nomials in /i*/^. A solution can be found if every coefficient of the /i*/^-terms is zero. Using 
conditions given in equations fl^ [T^ allows directly to find the set of solutions. 

Moulton (1963) describes the method and gives explicitly the expressions for Vp. and 
1pp. for i < 15. The integration constants with the full expressions are provided in the next 
equations 



= \^P^ 

= 0> 



Tpj 0, 



= ^Mp{5{l-M2f-9Mp + 5Ml), 



fpg — 0, 



r. 



r, 



9 

— cos 2r, 



PlO P 

2 T,j2 



rpu = 



^pi2 = ;^^p(175(l - - 535(1 - M2YM2 + 711(1 - M2yM^ - 535(1 - M2)M^ + 175M^) 
768 

- ^ ( 1 - 2M2) (3 cos r + - cos 3r) , 
2 9 



rp^3 = -Mp cos 2r, 



3 

32 " ■ 64 
175 



+ ^(25(1 - Ma)^ - 61A'/p + 25A'/|) cos 2r + 



^Q2A^p((^ - ^2)' - Mp + M|) cos4r, 
-^Mp(l - 2M2)(cosr + ^ cos3r), (44) 



and for the angle we have, 



V'P4 


= 0, 


^P5 


= 0, 


^P6 


= 0, 


^P7 


= 0, 
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= 0' 
3 

V^pio = -^Mp sin 2t, 
^Pii = 0, 

V^Pi2 = -|Mp(l-2M2)(sinr + ^sin3r), 
o y 

3 

^Pi3 = — Mpsin2r, 

5 35 
^Pi4 = — Mp((l-M2)' -4Mp + M|) sin 2r + —Mp((l-M2)' -Mp + M^) sin 4r, 

1 25 
= -Mp(l-2M2)(9sinr--sin3r), (45) 

where the integration constant for the radial equation [39] is calculated using the solution for 
six orders ahead. While for the angular equation HQ] the integration constant is obtained 
from the solution for ten orders ahead. 

It is worth to mention that the approximation to circumbinary periodic orbits given by 
Moulton is far from the solution as we have compared with the numerical solution, thus we 
have developed all the necessary terms in the expansion of the potential to reach a fairly 
good approximation to the numerical solution of the problem. 

In this way, an analytical expression is found for circumbinary orbits with high accuracy 
using the expressions HU HSl and the equations (HTl W2\ and [Tj). 



3.1. The Radial Extension of Circumbinary Discs (The Gap) 

Our purpose in this section is to provide a good estimation for the radius of the inner 
boundary of a circumbinary disc. Then we are looking for the closest stable orbits to the 
binary system. In this case, it is required to calculate orbits decreasing in size until con- 
secutive orbits intersect each other. The procedure we follow was to find a few new terms 
in the approximate solution for the disturbed orbit and calculate the larger inner orbit that 
intersects a contiguous one with the method described in Section [2^2] In this case, as in the 
circumstellar disc calculation (see Section [2?2]) . the angular correction ipp- is proportional to 
sinr, the intersections occur in the same manner at r = and r = tt, i.e., in the line joining 
the stars. The larger the number of terms considered for the approximation, the closer the 
orbit is to the orbit calculated in PSA05. However, the difference in sizes is large in spite 
of taking into account the terms up to /i^^^^ = (^)-2i/2 ^^j^g^^ gives high precision in the 



- 15 - 



circumstellar discs case. Table 2 gives for a set of values for Mi, the non-perturbed radius 
ttnp, the physical radius of the orbit in the intersection with the line connecting the stars, 
r(r = 0) and r(r = tt), and the analogous values in PSA05. Here, for each Mi there is an 
intersection at r = or at r = tt. Note that, for example, the system with Mi = 0.2 is the 
same system with Mi = 0.8, only rotated an angle tt, then, the intersection at a„p = 1.3920 
in r = in the former one, corresponds to a„p = 1.3920 in r = tt in the latter. Thus, Table 
2 gives radii for Mi < 0.5. 



3.2. Stability Analysis 

In the numerical approach, the solution for a closed orbit is searched by successive 
iterations. This means that an unstable orbit (surrounded by orbits far from the solution) is 
very hard to find. Unlike numerical calculations, analytically, one can calculate either stable 
and unstable orbits without any possible identification. One has to apply another criteria to 
look for the long-lived orbits as the ones in real discs, this criteria is taken from a stability 
analysis. 

In the case of circumstellar discs, the criteria to end up a disc was the intersection of 
orbits. In the circumbinary discs case, this is not applicable since in general orbits start 
being unstable before any intersection of orbits. 

Message (1959) defines an orbit as stable or unstable using the equation that describes 
normal displacements from a periodic orbit in the three-body restricted problem. These 
displacements are solved with terms proportional to exp(icr), where c is given by 

- (4 + A_i + Xiy + 2(Ai - A_i)c + A_iAi - = 0, (46) 

where 



Table 2: Size of the central gap in the circumbinary disc 



Ml 


^np,T=0 


r(T = 0) 


r(r = tt) 


r(r = 0)psA 


r(T = 7r)p5A 


0.1 


1.3547 


1.5660 


1.2956 


1.87 


1.80 


0.2 


1.3920 


1.6267 


1.3509 


2.04 


2.00 


0.3 


1.3733 


1.6179 


1.3637 


1.94 


1.90 


0.4 


1.3131 


1.5604 


1.3642 


1.92 


1.90 


0.5 


1.1643 


1.4183 


1.4183 


2.00 


2.00 
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A±i = Opo - 1 - {Po + (^±2)GpiGp-i - P±2P±3^1j^^I_i, (47) 

iy±i = -Ml^i, (48) 

and 9pg,9p^i are the main coefficients of the Fourier expansion of 6(r), where 6(r) is the 
function involved in the normal equation displacement, 



_1 + 0(r)g = 0, (50) 

which is given in Message (1959). The function 0(t) depends on the shape of the orbit, 
i.e. depends on r.lf for a specific orbit, c has an imaginary term then the orbit is unstable. 
Calculation of c is made for analytical orbits given by equations ([6] and [7j) with Vp and tpp 
estimated with equations (HTl H2|) and the coefficients given by equations HH and US Orbits 
smaller and larger that the intersecting-orbits given in Table 2 are considered. The values 
taken for the parameter a„p are given in Table 3. 

The parameter a = anp,psA corresponds to the analytical orbit closer to the innermost 
orbit of the circumbinary disc, found in PSA05, for several stellar mass values. Mi. Thus, if 
this set of orbits is (un)stable we expect that the numerical orbits also are (un)stable. We 
solve equation fH6|) for the orbits in Table 3 and the results are given in Table 4. There, the 
largest value of the imaginary part of c, looking at all the modes found, are given by Mi, 
ranging from 0.1 to 0.5. 

From Table 4, the trend for Max\Im{c)\ is to decrease when the orbit is moving away 
from the binary system, as expected. Note that Max\Im{c{anp,psA))\ is quite small, conse- 
quently, this value can be taken as a threshold for instability and the orbit in PSA05 will 



Table 3: Values of the parameter a^p for 



Ml 


(^np,l 


(^np, Inter. 


0.1 


1.25 


1.3547 


0.2 


1.29 


1.3920 


0.3 


1.27 


1.3733 


0.4 


1.21 


1.3131 


0.5 


1.06 


1.1643 



orbits studied in the instability analysis. 



(^np,2 


(^np,3 


Ctnp,PSA 




1.45 


1.65 


1.80 


2.75 


1.49 


1.69 


1.99 


2.79 


1.47 


1.67 


1.87 


2.77 


1.41 


1.61 


1.86 


2.71 


1.26 


1.46 


1.96 


2.56 
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naturally define the boundary between the unstable (not possible for gas orbits) and the 
stable (possible) part of the disc. 

It is worth noticing that this analysis can be applied either to numerical or analytical 
orbits. Only the Fourier coefficients of a function that depends on the orbit are required. 
Thus, the stability analysis described here is general, and can be used with any orbit of the 
restricted circular three-body problem. 



4. Test: Comparison with Numerical Results 



In this section we present a simple comparison between the results presented in this 
work and precise numerical results performed with the usual methods to integrate the motion 
equations in a circular binary potential. We show for this purpose two tests, the first is a 
direct qualitative comparison between numerical and the analytical approximation to some 
representative discs. The second including velocities, are the correspondent rotation curves. 



4.1. The Numerical Method 



The equations of motion of a test particle are solved for the circular restricted three 
body problem using an Adams integrator (from the NAG FORTRAN library). Cartesian 
coordinates are employed, with origin at the center of mass of the binary, in an inertial 
reference frame. Here we look for the families of circumstellar and circumbinary orbits at 
a given stellar phase of the stars (equivalent to look for the families in the non-inertial 
frame of reference). We calculate the Jacobi energy of test particles, in the non-inertial 
frame of reference co-moving with the stars, as a diagnostic for the quality of the numerical 
integration, conserved within one part in 10^ per binary period. 



Table 4: Maximum imaginary Part of c for the Orbits Studied in the Instability Analysis 

Max\Im{c)\\Mi 0.1 0.2 0.3 0.4 0.5 

c{anp,i) 2.65 X 10-3 5.64 x lO'^ 8.26 x 10"^ 9.38 x 10"^ 3.60 x 10"^^ 

c{anp,inters) 6.47 X 10"^ 1.43 X 10-^ 1.84 X 10-=^ 1.62 X 10-3 2.76 x 10-^^ 

c{anp,2) 2.35 X 10-^ 4.87 x 10-^ 5.83 x 10-^ 4.29 x 10-^ 4.35 x 10-^^ 

c(a„p,3) 4.96 X IQ-s 9.24 x IQ-^ 9.60 x 10"^ 5.53 x 10"^ 2.68 x IQ-^^ 

cianp,psA) 2.09 X 10-5 2.07 x 10-^ 2.72 x 10-^ 9.87 x 10-^ 

c{anp,A) 
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4.2. Geometrical Comparison of Discs 



In the case of binary systems is well known the shape and extension of discs are different 
from the single stars case, where the periodic orbits are circles with a radial unlimited (by 
the potential of the single star) extension. In binary systems the stellar companion exerts 
forces able to open gaps, limiting the extension of the discs and producing at the same time 
deformations in the shape of the discs making them visually different from single star cases. 
In this section we show a comparison between the discs obtained by the analytical approxi- 
mation provided in this work, the numerical solution and the keplerian approximation. 

To construct the analytical circumstellar discs we show in this section, we have employed 
the equations [6] and [7] that can be written in the inertial reference system. 



with r in the interval [0, 27i]. 

In the Figure [H we show the comparison for circumstellar discs between the approxima- 
tion presented in this work for periodic orbits in a binary system (left panels) and numerical 
calculations (right panels). The value for M2 is indicated on the top of the figure. The 
darker curve represents the Roche lobe. 

Likewise, in Figure [2] we show a qualitative comparison of a circumbinary disc between 
the analytical approximation (left panels) and numerical calculations (right panels). The 
distance of the stars with respect the center of mass of the binary is shown. 

The discs present a slight deformation of the orbits depending on the radius: the outer- 
most orbits are farer from being circles. The orbits are, a very good approximation, ellipses 
with an eccentricity depending on the radius to the central star, in the cincumstellar discs, 
and to the center of mass in the circumbinary discs case. Thus, we have calculated the m=2 
Fourier mode for all the orbits given by 




(51) 



(52) 





(53) 



i=0 



- 19 - 




-0.5 0.5 1 -0.5 0.5 1 



X (a) X (a) 

Fig. 1. — Circumstellar discs comparison between orbits with the approximation presented 
in this work and numerical results described in subsection 14.11 for the case M2 = 0.2. Left 
panel shows the analytical approximation provided in this work. Right panel shows the 
numerical approximation. The darker curve surrounding the discs is the Roche Lobe. The 
axes are in units of the distance between the stars, a 
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Ma =0.2 



M2=0.2 



-3 




- 



-4 



-4 



-2 





X (a) 



-2 





X id) 



Fig. 2. — Circumbinary discs comparison between orbits with the approximation presented 
in this work and numerical results described in subsection 14.11 for the case M2 = 0.2. Left 
panel shows the analytical approximation provided in this work. Right panel shows the 
numerical approximation. The axes are in units of the distance between the stars, a. The 
distance of the primary and secondary star with respect to the center of mass of the system 
is indicated by the black asterisks 
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where the index i refers to the number of evenly distributed (by interpolation) in angle points 
in a given loop, and is the total number of points in the loop, is the distance to 

the point i from the center of mass of the binary. In this manner, the average inner radius 
of a circumbinary disk is given in units of the semimajor axis a by the coefficient with 
k = 0. The ellipticity (ell) by ^Al + Bl with = 2 in the same units and transformed 
first to the ratio bi/ai, where a^, 6, are the semimajor and semiminor axis of the orbit i, and 
finally transformed to eccentricity (ecc) by. 



ecc = Vl - elP (54) 

which represents a more sensitive geometrical characteristic of the orbits than ellipticity. 
We have not considered in this comparison higher Fourier modes since they are negligible 
compared to the m=2 mode. 

In the Figure [3] we show the eccentricity vs the average radius for a) the orbits produced 
by the analytical approximation provided in this work (continuous lines in the 9 frames), 
b) the full numerical solution (open triangles), and for reference c) the keplerian (ecc = 0) 
solution (dashed lines). In this figure we show three different masses M2 = 0.01,0.2,0.4, 
at left, middle and right panels respectively (as indicated on the upper frames labels). The 
upper frames are referred to the circumprimary discs, the middle to the circumsecondary 
discs, and the lower are the circumbinary discs. The system of reference for each case is on 
the respective star for the circumstellar discs and in the center of mass for the circumbinary 
discs. 



4.3. Rotation Curves 

We also show here, for a quantitative comparison involving positions and velocities, the 
rotation curves for binary systems with M2 = 0.01, 0.2, 0.4. 

We have constructed the rotation curves by direct differentiation of equations [6] and [7] 
(or more precisely of equations [5T] and [521) . to obtain the total velocity {vc) at two different 
regions of the discs. For the primary discs we compared the rotation curves obtained along 
the X axis, on the left part of the discs (where the radial velocities are negative and prograde 
with the rotation of the binary). For the secondary we chose to compare the rotation curves 
along the x axis but on the right part of the disc (where the radial velocities are positive and 
prograde with the rotation of the binary). Once the velocities are calculated, we transformed 
the results to the inertial frame (in the circumstellar discs) by simply adding (or subtracting 
depending on the given disc) the velocity of the correspondent star, and adding or subtracting 



-22- 



M2 = 0.01 M2=0.2 M^ = OA 




0.2 0.4 0.6 0.1 0.2 0.3 0.4 0.1 0.2 0.3 




0.02 0.04 0.06 0.05 0.1 0.15 0.1 0.2 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 M 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


-, , , , 1 , , , , 1 , , , , 1 ,- 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



2 2.5 3 2 2.5 3 3.5 4 2 2.5 3 3.5 4 
R (Average Radius) 



Fig. 3. — Comparison between the orbital eccentricity vs radius produced by the analytical 
approximation provided in this work (continuous lines), the full numerical solution (open 
triangles), and for reference the keplerian (ecc = 0) solution (dashed straight lines). Three 
examples of different masses, M2 = 0.01, 0.2, 0.4, are indicated by the upper labels at left, 
middle and right panels respectively. The upper frames are referred to the circumprimary 
discs, the middle to the circumsecondary discs, and the lower are the circumbinary discs 
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also a factor fir, where fi is the angular velocity of the star and r is the radius of a given 
point in he rotation curve. The last because of the chosen reference system where this work 
solves the equations, which anchorages the discs to rotate with the system. 

In the Figure m and [5] we show the outer (top panels) and inner (bottom panels) regions of 
the primary and secondary rotation curve discs, respectively. By "outer" or "inner" regions 
we mean, in a reference system located in the primary (or in the secondary) star, measuring 
and angle from the line that joins the stars, the inner region corresponds to angle zero from 
this line, and the outer regions correspond to and angle of 180°. For both figures we plot 
the three different cases, M2 = 0.01,0.2,0.4. Three different types of lines indicate, a)the 
keplerian rotation curve on the top (filled circles); b) the analytical approximation provided 
in this work (continuous line); c) the numerical solution (open circles). The velocity and 
radius are given in code units (in such a way that G = 1, Mi + M2 = 1, a = 1 and fi = 1). 

The rotation curves are in approximately 70% of the total radius of the correspondent 
disc, keplerian curves. For the last 30% the keplerian discs velocity are systematically over 
the numerical result that solves with high precision the restricted three body problem. The 
analytical solution we provide here is very close to the numerical solution as expected. 

In the Figure [6] we show the rotation curves as in Figures H] and O but for the circumbi- 
nary disc. Although in this case both approximations (the one given in this work, and the 
numerical one) are very close to a keplerian rotation curve, it is not exactly keplerian. The 
velocities in the analytical approximation are sistematically under the keplerian curve and 
they are practically the same as the ones provided by the numerical solution, until the end 
(the beginning of the gap) where the numerical solution goes slightly over the analytical 
solution. We present only one case of mass ratio because other mass ratio would give very 
similar results. 

In the case of the full (numerical) solution, it is worth to notice we obtain in some cases 
discotinuities on the rotation curves and orbits mainly due to resonances. This is the case, 
for instance, in the Figure [6l where we appreciate a discontinuity very close to the resonance 
3:1 (r ~ —2.1a in the figure). The discontinuity is also noticeable in the Fig. [2] at the same 
radius. Due to the perturative fit performed in this work, finding this kind of discontinuities 
due to resonances, that result in general very narrow in radius but where the physics is 
highly non-linear is out of the reach of our approach. 
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5. CONCLUSIONS 

We have constructed an analytical set of equations based on perturbative analysis that 
result simple and precise to approximate the solution for the circular restricted three body 
problem. We choose some terms in the expanded equations of motion, which are relevant to 
the solution of the problem. The original equations are approximated with a set of equations 
that, in many cases, have an analytical solution. We show the goodness of our analytical 
approximation by direct comparison of geometry and rotation curves with the full numerical 
solution. 

The set of equations we prsent can provide any periodic orbit for a binary system, either 
circumstellar or circumbinary, and not only the outer edge radius of circumstellar discs or 
the inner edge radius (gap) of the circumbinary disc, without the need to solve the problem 
numerically. 

The relations provided are simple and straightforward in such a way that our approxi- 
mation can be used for any application where initial conditions of periodic orbits or complete 
periodic orbits are needed, or for direct study of the three body restricted circular problem. 
For instance, in order to introduce on a hydrodynamical or particle code, our initial condi- 
tions would result in much more stable discs than with keplerian initial conditions, or than 
by constructing the discs directly from hydrodynamical simulations of accretion to binaries 
that will require long times to obtain stable discs. Since our approximation is completely 
analytical it will have the obvious additional advantage of being computationally, extremely 
cheap, and easier to implement than the numerical solution. 

The periodic orbits respond uniquely to the binary potential and do not consider other 
physical factors such as gas pressure, or viscosity that work to build the fine details of 
discs structure. They are however, the backbone of any potential and their shape and 
behavior give the general discs phase space structure. In this manner, apart of all the 
possible hydrodynamical or particle discs applications, we can directly use them to study 
from the general discs geometry to rotation curves, or rarification and compression zones by 
orbital crowding, etc. 
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